#include "uhm/common.hxx"
#include "uhm/util.hxx"
#include "uhm/mat.hxx"

#include "uhm/mat/linal.hxx"

using namespace uhm;

#ifdef NOPIV 
typedef LDL_nopiv_Linal_  LDL_Linal;
#endif

#ifdef PIV 
typedef LDL_piv_Linal_    LDL_Linal;
#endif

#ifdef INCPIV 
typedef LDL_incpiv_Linal_ LDL_Linal;
#endif

int main (int argc, char **argv) {

  initialize(argc, argv);

  option_begin();
  option_end();

  Int_ mt = 5, mb = 3, rhs = 1, blk = 3;

  LDL_Linal A(Mat_Base_::DOUBLE, mt, mb, rhs, blk);
  A.create_flat_without_buffer();
  A.create_flat_buffer(true);
  A.create_hier();

  Double_ storage, max_storage;
  storage_counter(storage, max_storage);
  std::cout << "Storage = " << storage << ", Max used = " << max_storage << "  Bytes\n";

  random_spd_matrix(A);

  if (get_verbose())
    std::cout << A;

  A.decompose();

  copy(A.BT(), A.XT());
  A.XB().set_zero();

  A.solve_1_x();
  A.solve_2_x();

  if (get_verbose())
    std::cout << A;

  Double_ flop;
  A.decompose(flop);
  std::cout << "Decompose = " << flop << " FLOP\n";

  A.solve_1_x(flop);
  A.solve_2_x(flop);
  std::cout << "Solve = " << flop << " FLOP\n";

  A.free_all();
  
  finalize();
  return true;
}


